function [ outputSamples ] = WattersonFading( inputSamples, SamplingFrequencykHz, ...
    DopplerSpreadHz, pathDelaysMs, fadingMethod, Oversample, Fractional_Flag)
%WATTERSONFADING Summary of this function goes here
%   Detailed explanation goes here

numInputSamples=length(inputSamples);
outputSamples=zeros(1,numInputSamples);

%Make sure the path delays are sorted
pathDelaysMs=sort(pathDelaysMs);

%Make sure the first path delay is 0
pathDelaysMs=pathDelaysMs-pathDelaysMs(1);

[numRows, numColumns]=size(inputSamples);
if numRows >1
    % input is a column. make it row
    inputSamples=transpose(inputSamples);
end

pathDelaysSamples=round(pathDelaysMs*SamplingFrequencykHz);
numPaths=length(pathDelaysMs);

averageGain=zeros(1,numPaths);

% fadingMethod:
% 0 = No fading
% 1 = Uniform Power
% 2 = exponentilal Power profile

%=========================================================
% Select a sampling frequency in Hz about 100 times Doppler. call it virtual_fs
% 
% Generate a a fading waveform the number of samples is
% numInputSamples * virtual_fs/ (SamplingFrequencykHz*1000)

% interpolate this fading waveform for
% (SamplingFrequencykHz*1000)/virtual_fs in between virtual samples
%=========================================================


if fadingMethod==0 % No fading
    outputSamples=inputSamples;
    channelGain_out=ones(1,numInputSamples);
    return
elseif fadingMethod==1 % uniform Power Profile
    fadingTapPower=zeros(1 , 1+pathDelaysSamples(end) );
    fadingTapPowerDoppler=ones(1,numPaths);
    fadingTapPower(1+pathDelaysSamples)=fadingTapPowerDoppler;
elseif fadingMethod==2
    fadingTapPower=zeros(1 , 1+pathDelaysSamples(end) );
    fadingTapPowerDoppler=exp(-(0:numPaths-1)/numPaths);
    fadingTapPower(1+pathDelaysSamples)=fadingTapPowerDoppler;
else
    outputSamples=inputSamples;
    channelGain=1;
    return
end
fadingTapPower=fadingTapPower/sum(fadingTapPower);
fadingTapPowerDoppler=fadingTapPowerDoppler/sum(fadingTapPowerDoppler);
fadingTapGain=sqrt(fadingTapPower);
fadingTapGainDoppler=sqrt(fadingTapPowerDoppler);

%======================================================
if DopplerSpreadHz==0

    fadingGains=fadingTapGain .* (randn(1,length(fadingTapGain))+ ...
                                j*randn(1,length(fadingTapGain))) * sqrt(0.5) ;
    %fadingGains=1;
    % Make the fading symbol spaced
    % we can also make it more or less than one symbol, manually
    outputSamples= conv(inputSamples,fadingGains);
    channelGain=fadingGains;
    return
end
%======================================================

% reductionFactor=100;
% reducedFsHz=SamplingFrequencykHz*1000/ reductionFactor ;
% muDTm=DopplerSpreadHz/reducedFsHz ;

muDTm=0.01 ;
reducedFsHz=DopplerSpreadHz/ muDTm ;
reductionFactor=round(SamplingFrequencykHz*1000/reducedFsHz);

Lpath=(-100 : 100);
FilterTapGains=sqrt(sqrt(2 * pi * muDTm^2 )) * exp( -(pi*Lpath*muDTm).^2 );

filterLength=length(FilterTapGains);

% overSamplingFactor= round(SamplingFrequencykHz*1000/virtual_fs) ;

numFadingSamples=2*filterLength+ round(numInputSamples / reductionFactor)+ 1000; % just many samples

gaussianSamples=(randn(numPaths, numFadingSamples)+ j*randn(numPaths, numFadingSamples))*sqrt(0.5);

for pathCounter=1:numPaths
    tapGains=fadingTapGainDoppler(pathCounter) * conv(FilterTapGains , gaussianSamples(pathCounter,:) );
    
    %tmp=tapGains(2*filterLength:2*filterLength+round(numInputSamples/overSamplingFactor)+20);
    tmp=tapGains(2*filterLength:end);
    trueIndices=(1 : 1/reductionFactor : (100+numInputSamples/reductionFactor) );
    tapGains=interp1(tmp, trueIndices, 'linear', 'extrap');
    truncTapGains=tapGains(1:numInputSamples);
    
%     virtualIndex=overSamplingFactor*(0:length(tmp)-1)+1;
%     truncTapGains=interp1(virtualIndex, tmp, (1:numInputSamples),'linear','extrap');
    
%     tmp=interp(tmp,overSamplingFactor);
%     truncTapGains=tmp(1:numInputSamples);
    
    channelGain(pathCounter, 1: length(truncTapGains))=truncTapGains;
    
    %outputSamples=outputSamples+ truncTapGains .* circshift(inputSamples, [1, (pathCounter-1)*Oversample] );
    outputSamples=outputSamples+ truncTapGains .*  ...
                             [zeros(1,pathDelaysSamples(pathCounter)) , ...
                              inputSamples(1 : numInputSamples -pathDelaysSamples(pathCounter))];
end

if Fractional_Flag == 0
    numPaths_out = round(pathDelaysSamples(end)/Oversample);
else
    numPaths_out = round(pathDelaysSamples(end));
end

% channelGain_out(1,:) = channelGain(1,:);
% 
% channelGain_out(2:numPaths_out,:) = zeros(numPaths_out-1,length(channelGain));
%     
% channelGain_out(numPaths_out+1,:) = channelGain(end,:);


end








